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Abstract 

In this paper, we will analyze the behavior of the Neighbor Joining algorithm on five 
taxa and we will show that the partition of the sample (data) space for estimation of a tree 
topology with five taxa into subspaces, within each of which the Neighbor Joining algorithm 
returns the same tree topology. A key of our method to partition the sample space is the 
action of the symmetric group 6*5 on the set of distance matrices by changing the labels of 
leaves. The method described in this paper can be generalized to trees with more than five 
taxa. 



1 Introduction 



The Neighbor Joining (NJ) algorithm is introduced by ISaitou and Neil 19871 ] and it is widely 
used because of its accuracy and computational speed. A number of attempts have been made 
to understand the good results obtained with the Neighbor Joining algorithm, especially given 
the problems with the inference procedures used for estimating pairwise distances. For example, 
Brvantl 20051 ] showed that the Q-criterion (defined in ( [Ql ) in Section [2.2p is in fact the unique 
selection criterion such that it is linear, permutat ion equivariant, and cons istent, i.e. it correctly 
finds the tree corresponding to a tree metric and Gascuel and Steel 20061 ] showed a nice review 
of how the NJ algorithm works. 

One of the most important questions to study the behavior of the NJ algorithm is to analyze 
its performance with pairwise distances that are not tree metrics, especially when all pairwise 
distan ces are estimat ed via the maximum likelihood estimation (MLE). The following theorem 
due to lAttesonI [l999l ] is one of the results to address this question: 



Theorem 1.1. (Atteson (l999l ]) Neighbor- joining has Ic 



radius \ . 



This means that if the distance estimates are at most half the minimal edge length of the tree 
away from thei r true valu e theii the Neighbor Joining algorithm will reconstruct the correct 
tree. However, iLevv et al\ |2006l ] noted that Atteson's criterion frequently fail s to be satisfied 



even though the NJ algorithm r eturns the true tree to pology. Recent work of lMihaescu et al 



20061 ] extended Atteson's work. Mihaescu et al. 20061 ] showed that the NJ algorithm returns 



the true tree topology when it works locally for the quartets in the tree. This result gives another 
criterion when then NJ algorithm returns the correct tree topology and Attes on's theorem is 



a spe cial case of Mihaescu et al's theorem. The results of lAttesonl |1999l ] and iMihaescu et al 



20061 ]. however, depend on an input data and to analyze the performance of the NJ algorithm for 



any input data we need to study the sample space (the space of all possible distance matrices). 
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Therefore, this justifies to study the sample space which is the subspace of the vector space 
IR^a), where n is the number of taxa. In this paper, we will show the partition of the sample 
space for estimation of a tree topology with five taxa {0, 1, 2, 3, 4} into subspaces, within each of 
which the NJ algorithm returns the same tree topology. For five taxa, the sample space is the 
subspace IR^ of the vector space IRelR'^elR^ Note that our method to partition the sample 
space into subspaces described in this paper can be generalized to trees with more than five 
taxa. 

For every input distance matrix, the NJ algorithm returns a certain tree topology. It may 
happen that the minimum Q-criterion is taken by more than one pair of taxa at some step. 
In practice, the NJ algorithm will then have to choose one cherry in order to return a definite 
result, but for our analysis we assume that in those cases the NJ algorithm will return a set 
containing all tree topologies resulting from picking optimal cherries. 

There are only finitely many tree topologies, and for every topology t we get a subset Dt of 
the sample space (input space) such that for all distance matrices in Dt one possible answer of 
the NJ algorithm is t. We aim at describing these sets Dt and the relation between them. 

It turns out that we can break up each Dt into easier subsets, by considering not only the 
topology returned by the NJ algorithm but also the order in which cherries are picked. Consider 
the tree in figure [T^a). There are two possible ways in which the NJ algorithm may reconstruct 
this tree, namely by first picking the cherry 1-0 or by first picking the cherry 4-3. The set of 
input distance matrices (vectors in IR© IR^ © IR^) for which the NJ algorithm will pick exactly 
the same cherries in the same order then it turns out to be a cone, so it can easily be described 
by a set of halfspaces. 

One important tool in our analysis will be the action of the symmetric group S5 on the set 
of distance matrices by changing the labels of leaves. Because there is only one unlabeled tree 
topology in the case of five taxa, any two labeled tree topologies can be mapped onto each other 
by changing leaf labels. Thus it suffices to analyze just one Dt- 



2 The Neighbor-Joining algorithm 
2.1 Input data 

The NJ algorithm is a distance based method which takes a distance matrix, a symmetric matrix 
{dij)o<ij<n-i with da = representing pairwise distances of a set of n taxa {0, 1, ... ,n — 1}, 
as the input. Through this paper, we do not assume anything on an input data except it is 
symmetry and da = 0. Because of symmetry, the input can be seen as a vector of dimension 
m := (2) = ^n(n — 1). We arrange the entries row-wise, so in the case n = 4 we get: 



( ° 


do 


di 


ds \ 
di 


do 





d2 


di 


d2 





d5 


\ ds 


di 


d5 


/ 



We denote row/column-indices by pairs of letters such as a, b, c, d, while denoting single indices 
into the "flattened" vector by letters i,j,.... The two indexing methods are used simultaneously 
in the hope that no confusion will arise. Thus, in the four taxa example we have do,i = di^ = do- 
In general, we get di = dab = dba with 



max 



k I -k{k -1) <i 



1 



+ 



1 



+ 2i 



-{a — l)a, 
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and for c > d we get 



2.2 The Q-Criterion 



d. 



■c,d — dc[c-l)/2+d- 



The NJ algorithm starts by computing the so called Q- criterion or the cherry picking criterion, 
given by the formula 



n-1 



n-1 



i-k,b- 



(Q) 



Qafi ■■= {n - 2)dafi - ^ da^k - '^dk 

This is a key of the NJ algorithm to choose which pair of taxa is a neighbor. 

Theorem 2.1. (ISaitou and Nei (l987l] : IStudier and Kepplerl Let da b for all pair 

of taxa {a, b} be the tree metric corresponding to the tree T. Then the pair {x.y} which minimizes 
Qa,b for all pair of taxa {a, b} forms a neighbor. 

The resulting matrix is again symmetric, and ignoring the diagonal entries we can see it as 
a vector of dimension m just like the input data. Moreover, the Q-criterion is obtained from the 
input data by a linear transformation: 

q = A(")d, 

and the entries of the matrix j4^'"^ are given by 



A, 



(n) 



A 



in) 
ab,cd 



n — 4 if i = j, 

-1 if ? / j and {a, b} n {c, d} / 0, 
else, 





-1 


-1 


-1 


-1 





\ 




-1 


-1 


-1 





-1 






-1 -1 








-1 


-1 






-1 -1 








-1 


-1 






-1 


-1 


-1 





-1 




v 


-1 


-1 


-1 


-1 





/ 



where a > 6 is the row/column- index equivalent to i and likewise for c > d and j. When no 
confusion arises about the number of taxa, we abbreviate A^""^ to A. In the case of four taxa, 
we get 



^(4) 



After computing the Q-criterion q, the NJ algorithm proceeds by finding the minimum entry 
of q, or, equivalently, the maximum entry of — q. The two nodes forming the chosen pair (there 
may be several pairs with minimal Q-criterion) are then joined ("cherry picking"), i.e. they 
are removed from the set of nodes and a new node is created. Suppose out of our Ti taxa 
{0, . . . ,n — 1}, the first cherry to be picked is m — 1, so the taxa n — 2 and n — 1 are joined 
to form a new node, which we view as the new node number n — 2. The reduced pairwise 
distance matrix is one row and one column shorter than the original one, and by our choice of 
which cherry we picked, only the entries in the rightmost column and bottom row differ from 
the original ones. Explicitly, 



d'i 



di 



for < i < 



\{di + (ij+(„_2) - dm-i) for ("2 ) <i < 



(n~l\ 
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and we see that the reduced distance matrix depends hnearly on the original one: 

d' = Rd, 

with R = (rij) E ir{'"-"+i)x'", where 

for < i = j < 



for (-2) < ^ < 



1 

1/2 

1/2 for ("-2) < i < 
-1/2 for ("-^ < i < 
else 



\ 2 



i + n-2 
m — 1 



The process of picking cherries is repeated until there are only three taxa left, which are then 
joined to a single new node. 



3 The cones Cab,c 

In the case of five taxa there is just one unlabeled tree topology (cf. figured]) and there are 15 
distinct labeled trees: We have five choices for the leaf which is not part of a cherry and then 
three choices how to group the remaining four leaves into two pairs. For each of these labeled 
topologies, there are two ways in which they might be reconstructed by the NJ algorithm: There 
are two pairs, any one of which might be chosen in the first step of the NJ algorithm. 

2 




Figure 1: (a) A tree with five taxa (b) The same tree with all edges adjacent to leaves reduced 
to length zero. The remaining two edges have lengths a and /?. 

For distinct leaf labels a, b and c G {0,1,2,3,4} we define Cab,c to be the set of all input 
vectors for which the cherry a-b is picked in the first step and c remains as single node not part 
of a cherry after the second step. For example, the tree in figure [Dj^a) is the result for all vectors 
in Cio,2 U C43^2- Since for each tree topology ((a, b),c, {d, e)) (this tree topology is written in the 
Newick format) for distinct taxa a,b,c,d,e € {0,1,2,3,4}, the NJ algorithm returns the same 
tree topology with any vector in the union of two cones Cab,c U Cde,c, there are 30 such cones in 
total, and we call the set of these cones C. 



3.1 The shifting lemma 



We first note that there is an n- dimensional linear sub space of IR"* which does not affect the 
outcome of the NJ algorithm (see Mihaescu et al . 2006l |). For a node a we define its shift vector 
Sa by 

Jl if a G {b,c} 
I else 



(Sa)fe,c 
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which represents a tree where the leaf a has distance 1 from ah other leaves and all other distances 
are zero. The Q-criterion of any such vector is —2 for all pairs, so adding any linear combination 
of shift vectors to an input vector does not change the relative values of the Q-criteria. Also, 
regardless of which pair of nodes we join, the reduced distance matrix of a shift vector is again 
a shift vector (of lower dimension), whose Q-criterion will also be constant. Thus, for any input 
vector d, the behavior of the NJ algorithm on d will be the same as on d + s if s is any linear 
combination of shift vectors. We call the subspace generated by shift vectors S. 

We note that the difference of any two shift vectors is in the kernel of A, and the sum of all 
shift vectors is the constant vector with all entries equal to n. If we fix a node a then the set 

{sa - Sft I 6 / a} 

is linearly independent. 



3.2 Permuting leaf labels 

Because there is only one unlabeled tree topology, we can map any labeled topology to any other 
labeled topology by only changing the labels of the leafs. Such a change of labels also permutes 
the entries of the distance matrix: For example, the result of swapping leaves and 1 is the 
distance vector 

d' := {do, d2, di,d4, ds, 4, ^7, de, ds, dg)'^ . 

In this way, we get an action of the symmetric group 55 on the input space R^^, and the 
permutation a £ S5 maps the cone Cab.c linearly to the cone Cf^(^a)a{b),cr{c)- Therefore any 
property of the cone Cab,c which is preserved by linear transformations must be the same for all 
cones in C, and it suffices to determine it for just one cone. 

The action of S5 on R^^ decomposes into irreducible representations by 

]R©]R^© K^, 
=s =.w 

where the first summand is the subspace of all constant vectors and the second one is the kernel 
of The sum of these two subspaces is exactly the space S generated by the shift vectors. 

The third summand, which we call W , is the orthogonal complement of S and it is spanned by 
vectors Wab.cd in W with 



{Wab,cd) 



xy ■- 



1 if xy = ab or xy = cd 
— 1 xy = ac or xy = bd 
else 



where a, b, c and d are pairwise distinct taxa in {0, 1, 2, 3, 4} and {wab,cd)xy is the x-yt\\ coordi- 
nate of the vector Wab,cd- One linearly independent subset of this is 



Wl 



■^34,12, 



■^40,23- 



Note that the 5-cycle (01234) of leaf labels cyclically permutes these basis vectors, whereas the 
transposition (01) acts via the matrix 



r := - 
2 



/ 2 





Vo 



1 \ 

-1 
-1 
-1 

1 / 
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Because a five-cycle and a transposition generate S^, in principle this gives us complete infor- 
mation about the operation. 

3.3 The cone C43 2 

Since we can apply a permutation a £ S^j, without loss of generality, we suppose that the first 
cherry to be picked is the cherry 9, which is the cherry with leaves 3 and 4. This is true for all 
input vectors d which satisfy 

(h9,i,d) >0fori = 0,...,8, 

where the vector 

hff := -AW(e.-e,) 

is perpendicular to the hyperplane of input vector for which cherries i and j have the same 
Q-criterion, pointing into the direction of vectors for which the Q-criterion of cherry i is lower. 

We let ri, r2 and r3 be the first three rows of —A^^^R^^\ If (ri,d) is maximal then the 
second cherry to be picked is 0-1, leaving 2 as the non-cherry node, and similarly r2 and lead 
to non-cherry nodes 1 and 0. This allows us to define the set of all input vectors d for which 
the first picked cherry is 3-4 and the second one is 0-1: 

Cu,2 := {d I (hg,,, d) > for i = 0, . . . , 8, and (rg - ri, d) > 0, (rg - r2, d) > 0}. (2) 

We have defined this set by 11 bounding hyperp lanes. However, in fact, the r esulting cone has 



only nine facets. A computation using polymake [Gawrilow and Joswigi . |2000| | reveals that the 
two hyperplanes hg 1 and 119^2 are no longer faces of the cone, while the other nine hyperplanes 
in ([2]) give exactly the facets of the cone. 

3.4 The Rays of C 



Again using polymake [Gawrilow and Joswigj . l2000l | we find that 6*43^2 is the positive cone 



spanned by fourteen rays. The union of the orbits of these rays under the action of ^5 is a 
set of 82 vertices, which we call TZ. Each of the cones in C is spanned by a certain subset with 
14 vertices each of TZ. 

The set TZ has three orbits under the action of the symmetric group 6*5 . We characterize the 
rays in these orbits using the graphs in Figure [2j In these graphs, nodes which are connected 
by an edge form a pair with the minimum Q-criterion in the first step of the NJ algorithm, and 
the labels to each edge show which nodes are possible as the remaining unpaired node. 

b ^^'i 



(a) b^c (b) y (c) 

/ \ a I be be 



cd/ I \bc 
bd \ 



bd \ ^ ^ cd \ I / be 

e ^ e 



Figure 2: Graphs characterizing the rays in TZ. Here a, b, c, d and e are pairwise distinct labels 
between and 4. 

For each graph in figure [2] and each assignment of the five leaves 0, . . . , 4 to the variables 
a, . . . , e we get a set G of vectors which belong exactly to the cones indicated by the graph. By 
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what we have said in sections 13.11 and 13.2^ the set G has the form 

G = {ag + s I a > and s G 5}, 

where S is the set of all shifting vectors in the sample space (which is the subspace in R^^). 
Thus C is described by a vector g which is unique up to normalization by a positive constant. 
We call this vector g^^),^g, Sabcde siftcde' depending on which graph it corresponds to. 
Some of these graphs coincide for different choices of a, . . . , e. 

• The graphs in EJa) are strings of four taxa. There are five ways to choose the four taxa, 
and f = 12 ways in which to arrange them (reversing the order of the string does not 
matter). In total there are 60 rays of this type. 

fa) 

Thus for example the string 0-1-2-3 refers to an input vector goi234 which lies in the 
intersection of the five cones Cio,4, C2i,0i (^21,35 (^21,41 ^30,4 but not in any other cone. It 
is worth noting that this vector is a ray of the cones 6*21,0 > 6*21,3 and 621,4 but not of the 
cones Cio,4 and 632,4, which is indicated by the dotted edges in figure [2j This implies that 
the cones in C do not form a fan, and in particular no polytope in IR^^ can have C as its 
normal fan. 

• There are ^ possibilities for the cycle graph in figure [2]J^b) , each of which represents an 
input vector which forms a ray of ten of the cones in C. 

• Finally, there are (2) = 10 graphs of type[2|^c), because these are determined by the two 
nodes we choose for a and e. Each of these corresponds to an input vector that forms a 
ray of twelve cones. 

If we collect all those rays in TZ which are rays of 6*43,2 we get 

• 6 rays of type (a): 

(a) (a) (a) (a) (a) (a) 

§03412 ' §03421 ' §13402 ' §13420 ' §23401 ' §23410 



4 of type (b): 
4 of type (c): 



(b) (b) (b) (b) 

§34201 ' §34210 ' §34012 ' §34102 



(c) (c) (c) (c) 

§34201 ' §34210 ' §34021 ' §34120 • 



We call these vectors gi, . . . , gi4. Any vector v G 643,2 is a linear combination of these fourteen 
vectors with nonnegative linear coefficients, say 



V = ^ aigi 



where we number the fourteen rays arbitrarily. By linearity, in both steps of the NJ algorithm, 
the Q-criteria are linear combinations of the Q-criteria for the gi, using the same linear co- 
efficients ai. Furthermore, because all of the vectors gj lie in 643,2, in both steps of the NJ 
algorithm there is a pair of leaves such that the minimum of the Q-criteria is attained at this 
pair for all gj, and therefore the minimum Q-criterion of the linear combination ^ ajgj is the 
linear combination of the minimum Q-criteria of the individual gj. This minimum is attained 
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exactly by those pairs which have minimum Q-criterion in all gj with strictly positive linear 
coefficient a^, we get: 

V G Cab,c iff Si G Cab,c for all i with > 0. 

Therefore the rays in TZ are the input vectors for which the NJ algorithm is least stable. We 
can give explicit descriptions of these vectors using the graphs in Figure [3l The graphs give the 
distance value assigned to each pair of leaves. For example, we get 




Figure 3: Graphs describing the vectors gabcde for each of the three types of graphs. Labels 
assigned to edges denote the distance assigned to the corresponding pair of leaves in the vector, 
and edges of the same type receive the same label. 

3.5 The Cone of Tree Metrics 

Suppose our input vector is a tree metric. By relabeling the tree and using the shifting lemma, 
we can reduce the tree to Figure [T]^b), where all edges connecting to leaves have length zero. 
The edge-lengths a and /3 determine the distance vector v as follows: 

V = a(0, 1, 1, 1, 1,0, 1, 1,0, 0)"^ + /?(0,0,0, 1, 1, 1, 1, 1, 1,0)^ =: aa + /3b. 

If we assume these to be non-negative, we get a cone spanned by the two vectors a and b. These 
satisfy 

a e Cio,2 n Cio,3 n Cio,4, b e 043,0 n 043,1 n C43,2 and ^(a + b) e Cio,2n 043,2, 

and none of these vectors lies in any other cone. 

We normalize to a + /? = 1, and because of symmetry we may assume that a > (3, or 
equivalently a > ^. Then we can write v as 

V = aa + (1 — a)h 

= - a+ (1 - a)b' 

Because both a and b' are in the convex set Cio,2 5 the convex linear combination v is also in 
Cio,2; proving correctness of the NJ algorithm for five taxa. 

We can also compute the I2 distance of v from the faces of the cone Cio,2- There are 9 
hyperplanes defining Cio,2; but we may ignore one of them, which defines the common facet 
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with 6*43^2 • For any hyperplane H such that a and b' he on the same side of H, the distance 
d{v, H) between H and v is given by 



(i(v, H) = (^a- d{si, H) + {1- a)d{h', H), 



and taking the minimum of this over the eight remaining faces of Cio,2 we obtain 

d(v, (Cio,2 n C743,2) = ) = —j=-, 

where (Cio,2n (743,2)'^ is the complement of Cio,2n 6*43,2 • If we divide this by the length f3 = 1 — a 
of the smaller of the two interior edges, we get an ^2-radius of l/Vs ^ 0.577. Note that because 
our method relies on orthogonal projections, we get £2 bounds instead of £00 bounds. 



4 Simulation results 

In this section we will analyze how the tree metric for a tree and pairwise distances estimated 
via the maximum likelihood estimation locate in the partition of the s ample space . Part icularly, 
we analyze subtrees of the two parameter family of trees described by Ota and Li . 2OO0I ]. These 
are trees for which the NJ algorithm has difficulty in resolving the correct topology. In order 
to understand how they locate to each other, we simulated 10,000 data sets on each of the two 
tree shapes, Ti and T2 (Figured]) at the edge length ratio, a/b = 0.0 3/0.42 for sequences of 
length 500BP under the Jukes-Cantor mo del Jukes and Cantor . 1969l |. We also repeated the 
ru ns with the Kimur a 2-parameter model Kimura . 1980t . Thev are the cases (on eight taxa) 



m 



^Ota and Li l2000l| that the NJ algorithm had most difficulties in their simulation study (also 

the same as in iLevv et all l2006l ]l. Each set of 5 sequences are generated via evolver from 



PAML package [Yang) . Il997l | under the given model, evolver generates a set of sequences given 
the model and tree topology using the birth-and-death process. For each set of 5 sequences, 
we compute first p airwise distances via the heuristic MLE method using a software fastDNAml 



Qlsen et all Il994i |. To compute cones, we used MAPLE and polymake. 



Figure 4: Ti and T2 tree models which are subtrees of the tree models in |Ota and Li l200nl | 
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To study how far each set of pairwise distances estimated via the maximum hkehhood esti- 
mation (which is a vector y in ]R^) locates from the cone, where the additive tree metric hes, in 
the sample space, we calculated the ^2-distance between the cone and a vector y. 

Suppose we have a cone C defined by hyperplanes ni, . . . , n^, i.e. 

C = {x I (rij, x) > for i = 1, . . . , r}, 

and we want to find a closest point in C from some given point v. Because C is convex, for 
^2-norm there is only one such point, which we call u. If v G C then u = v and we are done. If 
not, there is at least one rij with (nj,v) < 0, and u must satisfy (nj,u) = 0. 

Now the problem reduces to a lower dimensional problem of the same kind: We project v 
orthogonally into the hyperplane H defined by (nj,x) = and call the new vector v. Also, 
C n H is a facet of C, and in particular a cone, so proceed by finding the closest point in this 
cone from v. 

We say an input vector (distance matrix) is correctly classified if the vector locates in one 
of the cones where the vector representation of the tree metric (noiseless input) lies. We say an 
input vector is incorrectly classified if the vector locates in the complement of the cones where 
the vector representation of the tree metric lies. For input vectors (distance matrices) which 
are correctly classified by the NJ algorithm, we compute the minimum distance to any cone 
giving a different tree topology. This distance gives a measure of robustness or confidence in 
the result, with bigger distances meaning greater reliability. The results are plotted in the left 
half of Figure [5] and in Figure [6l Note that the distance of the noiseless input, i.e. the tree 
metric from the tree we used for generating the data samples, gives an indication of what order 
of magnitude to expect with these values. 
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Figure 5: Distances of correctly (left) and incorrectly (right) classified input vectors from the 
closest incorrectly/correctly classified vector. 
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Figure 6: Mean and variance of the distances of correctly classified vectors from the nearest 
misclassified vector. 
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For input vectors to which the NJ algorithm answers with a tree topology different from the 
correct tree topology, we compute the distances to the two cones for which the correct answer 
is given and take the minimum of the two. The bigger this distance is, the further we are off. 
The results are shown in the right half of Figure [5] and in Figure [71 
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Figure 7: Mean and variance of the distances of misclassified vectors to the nearest correctly 
classified vector. 



From our results in Figure [6] and Figure [71 one notices that the NJ algorithm returns the 
correct tree more often with T2 than with Ti. T hese results are consistent with the results in 
Gascuel and Steel . 20061 : Mihaescu et ai . 2006| . Note that any pos sible quartet in T has a 



smaller (or equal) length of its internal e dge than in T2 (see Fig ure [3]) . iGascuel and Steell 2006l | 



defined this measure as neighborliness. [Mihaescu et al\ |2006l | showed that the NJ algorithm 



returns the correct tree if it works correctly locally for the quartets in the tree. The neighborliness 
of a quartet is one of the most important factors to reconstruct the quartet correctly, i.e., the 
shorter it is the more difficult the NJ algorithm returns the correct quartet. Also Figure [5] shows 
that most of the input vectors lie around the border lines of cones, including the noiseless input 
vector (the tree metric). This shows that the tree models Ti and T2 are difficult for the NJ 
algorithms to reconstruct the correct trees. All source codes for these simulations described in 
this paper will be available at authors' websites. 
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